This tutorial is part of viskex.

import os
import urllib

import dolfinx.io
import dolfinx.mesh
import mpi4py.MPI
import numpy as np
import numpy.typing

import viskex

Read in and plot mesh

msh_filename = "data/garda.msh"
if not os.path.isfile(msh_filename):
    os.makedirs("data", exist_ok=True)
    msh_url = (
        "https://raw.githubusercontent.com/FEMlium/FEMlium/main/"
        "tutorials/01_introduction/data/garda.msh")
    with urllib.request.urlopen(msh_url) as response, open(msh_filename, "wb") as msh_file:
        msh_file.write(response.read())

mesh, subdomains, boundaries = dolfinx.io.gmshio.read_from_msh(
    "data/garda.msh", comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)

Info    : Reading 'data/garda.msh'...
Info    : 2390 entities
Info    : 3937 nodes
Info    : 7872 elements
Info    : Done reading 'data/garda.msh'

viskex.dolfinx.plot_mesh(mesh)

viskex.dolfinx.plot_mesh(mesh, dim=2)

viskex.dolfinx.plot_mesh(mesh, dim=1)

viskex.dolfinx.plot_mesh(mesh, dim=0)

Plot subdomains

viskex.dolfinx.plot_mesh_tags(mesh, subdomains, "subdomains")

Plot boundaries

viskex.dolfinx.plot_mesh_tags(mesh, boundaries, "boundaries")

viskex.dolfinx.plot_mesh_entities(mesh, 1, "boundaries", boundaries.indices[boundaries.values >= 1])

Plot a scalar field

scalar_function_space = dolfinx.fem.FunctionSpace(mesh, ("CG", 2))

centroid = np.mean(mesh.geometry.x[:, :2], axis=0)

def scalar_field_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.float64]:
    """Evaluate the scalar field."""
    rho = np.sqrt((x[0] - centroid[0])**2 + (x[1] - centroid[1])**2)
    theta = np.arctan2(x[1] - centroid[1], x[0] - centroid[0])
    return rho / np.sqrt(1 - 0.5 * np.cos(theta)**2)  # type: ignore[no-any-return]

scalar_field = dolfinx.fem.Function(scalar_function_space)
scalar_field.interpolate(scalar_field_eval)

viskex.dolfinx.plot_scalar_field(scalar_field, "scalar")

Plot a vector field

vector_function_space = dolfinx.fem.VectorFunctionSpace(mesh, ("CG", 2))

def vector_field_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.float64]:
    """Evaluate the vector field."""
    rho = np.sqrt((x[0] - centroid[0])**2 + (x[1] - centroid[1])**2)
    theta = np.arctan2(x[1] - centroid[1], x[0] - centroid[0])
    values = np.zeros((2, x.shape[1]))
    values[0] = - rho * np.sin(theta)
    values[1] = rho * np.cos(theta)
    return values

vector_field = dolfinx.fem.Function(vector_function_space)
vector_field.interpolate(vector_field_eval)

viskex.dolfinx.plot_vector_field(vector_field, "vector")

viskex.dolfinx.plot_vector_field(vector_field, "vector", glyph_factor=0.1)